   function [stt,etamat,smooth_st,yfor,vstar,logLnc,countmat,countobs,countszero,outsol]=lmj_kfilterall(funcmod,fullp,y,flag_mode,ins,varargin)
 % function [stt,etamat,smooth_st,yfor,vstar,logLnc,countmat,countobs,countszero,outsol]=lmj_kfilterall(funcmod,fullp,y,flag_mode,ins)
% =========================================================================
% KFILTER_ALL.M 
% 
% function [stt,etamat,smooth_st,yfor,vstar,logLnc,countmat,countobs,countszero]=kfilter_all(funcmod,fullp,y,flag_mode,ins)
%
% Input 
% FUNCMOD 
% FULLP 
% Y 
% FLAG_MODE  == 1 Shut all but onw 
%            == 2 One state at a time 
%            == [] No counterfactual 
% OUTSOL     Spits out solution matrices GG,RR,SDX,ZZ,CC in a structure 
% INS       Structure for filter initialization 
% 
% Last Modified March 6 2008 
% =========================================================================
if nargin < 5 || isempty(ins)
    flag_ins=0; 
    disp('No restrictions on initial state') 
else 
    flag_ins=1; 
    disp('Use INS for initial states'); 
end 
if nargin < 4 || isempty(flag_mode) || flag_mode==0 
    flag_counter=0; 
else
    flag_counter=1; 
    if flag_mode~=1 && flag_mode~=2 
        error('FLAG mode must be 1 or 2') 
    end 
    dispaj('Mode for the counterfactuals: ',flag_mode); 
end 
% [TT,CC,RR,eu,SDX,ZZ]=feval(funcmod,fullp);
[TT,RR,CC,eu,SDX,ZZ]=feval(funcmod,fullp, varargin{:});
outsol.GG=TT; 
outsol.RR=RR; 
outsol.CC=CC; 
outsol.eu=eu; 
outsol.SDX=SDX; 
outsol.ZZ =ZZ; 


if ~isequal(eu,[1;1]);error('Model is not determinate');end
if any(CC~=0)
    disp('Model demeaning the data');
    ydem=y-repmat((ZZ*CC)',[size(y,1) 1]);
else
    ydem=y;
end
switch flag_ins
    case 0
        [stt,etamat,smooth_st,yfor,vstar,logLnc,countszero]=kfilter_full(ydem,ZZ,TT,RR,SDX);
    case 1
        [stt,etamat,smooth_st,yfor,vstar,logLnc,countszero]=kfilter_full(ydem,ZZ,TT,RR,SDX,ins);
end
stt=stt';
if flag_counter==0
    countmat=[];countobs=[];
    return
end
nobs=size(y,1);
nx=size(SDX,1);
nst=size(smooth_st,2);

disp('Begin Counterfactual')
countmat=zeros(nobs,nst,nx);
countobs=zeros(nobs,size(ZZ,1),nx);
% Begin Counterfactual
for ii=1:nx
    tempeta=zeros(nx,nobs);
    if flag_mode==2
        tempind=(1:nx);tempind(ii)=[];
        tempeta(tempind,:)=etamat(:,tempind)';
    else
        tempeta(ii,:)=etamat(:,ii)';
    end
    % ==================================================
    % Loop follows KFILTER_FULL.m
    ysim=zeros(size(ZZ,1),nobs);
    stsim=zeros(nst,nobs);
    stsim(:,1)=countszero(:,ii);
    ysim(:,1)=ZZ*stsim(:,1);
    for jj=2:nobs;
        stsim(:,jj)=TT*stsim(:,jj-1) + RR*(tempeta(:,jj));
        ysim(:,jj)=ZZ*stsim(:,jj);
    end
    % =====================================================
    countmat(:,:,ii)=stsim(:,:)';
    countobs(:,:,ii)=ysim';
end
% Check that the sum of the counterfactual matrices adds up to the
% smoothed states October 23 2007
% If FLAG_MODE == 2 the counterfactual need not recover the observables
if flag_mode==1
    checkysim=max(abs(sum(countobs,3)-ydem));
    disp('FLAG_MODE==1 Check Counterfactuals')
    dispaj('Max difference counterfactual & observables is:',max(checkysim));
    if max(checkysim) > 1e-4;
        error('Cannot recover original states with counterfactual');
    end
    checksim=max(abs(sum(countmat,3)-smooth_st(:,:)));
    dispaj('Max difference counterfactual & smoothed states:',max(checksim));
    if max(checksim) > 1e-4;
        error('Cannot recover original states with counterfactual');
    end
else
    disp('FLAG_MODE==2 hence will not recover observables')
end 